Here are the details of the Bayesian multi-level modelling. Our general approach is:
Define Priors
In this section we will specify some priors. We then then use a prior-predictive check to assess whether the prior is reasonable or not (i.e., on the same order of magnitude as our measurements).
Fixed Effects
Our independent variable is a categorical factor with 16 levels. We will drop the intercept from our model and instead fit a coefficent for each factor level (\(y \sim x - 0\)). As our dependant variable will be log-transformed, we can use the priors below:
prior <- c(
set_prior("normal(0,2)", class = "b"),
set_prior("cauchy(0,2)", class = "sigma"))
Group-level Effects
We will keep the default weakly informative priors for the group-level (‘random’) effects. From the brms documentation:
[…] restricted to be non-negative and, by default, have a half student-t prior with 3 degrees of freedom and a scale parameter that depends on the standard deviation of the response after applying the link function. Minimally, the scale parameter is 10. This prior is used (a) to be only very weakly informative in order to influence results as few as possible, while (b) providing at least some regularization to considerably improve convergence and sampling efficiency.
Prior Predictive Check
Now we can specify our Bayesian multi-level model and priors. Note that as we are using sample_prior = 'only', the model will not learn anything from our data.
m_prior <- brm(data = d_eeg,
rms ~ wg-1 + (1|subject),
family = "lognormal",
prior = prior,
iter = n_iter ,
sample_prior = 'only')
We can use this model to generate data.
## Observations: 320,000
## Variables: 2
## $ key <chr> "P2", "P2", "P2", "P2", "P2", "P2", "P2", "P2", "P2", "P2"…
## $ value <dbl> 1.98845697, 9.26489539, 9.71836157, 0.07553109, 0.27777647…
We can see that i) our priors are relatively weak as the predictions span several orders of magnitide, and ii) our empirical data falls within this range.
Compute Posterior
Fit Model to EEG Data
We will now fit the model to the data.
## Family: lognormal
## Links: mu = identity; sigma = identity
## Formula: rms ~ wg - 1 + (1 | subject)
## Data: d_eeg (Number of observations: 400)
## Samples: 4 chains, each with iter = 10000; warmup = 5000; thin = 1;
## total post-warmup samples = 20000
##
## Group-Level Effects:
## ~subject (Number of levels: 25)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.38 0.06 0.28 0.51 1.01 1426 1756
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## wgCM -0.48 0.09 -0.66 -0.31 1.00 869 1662
## wgCMM -0.12 0.09 -0.30 0.06 1.00 928 1848
## wgP2 -0.95 0.09 -1.12 -0.77 1.00 916 1827
## wgP3 -0.82 0.09 -0.99 -0.64 1.00 928 1695
## wgP31M -0.43 0.09 -0.61 -0.25 1.00 893 1772
## wgP3M1 -0.14 0.09 -0.31 0.04 1.00 870 1695
## wgP4 -0.48 0.09 -0.66 -0.31 1.00 921 1796
## wgP4G -0.25 0.09 -0.43 -0.08 1.00 911 1812
## wgP4M 0.15 0.09 -0.02 0.33 1.00 894 1818
## wgP6 -0.48 0.09 -0.66 -0.30 1.00 891 1783
## wgP6M -0.04 0.09 -0.22 0.13 1.00 866 1708
## wgPG -0.93 0.09 -1.11 -0.76 1.00 873 1713
## wgPGG -0.72 0.09 -0.90 -0.55 1.00 886 1660
## wgPM -0.59 0.09 -0.77 -0.42 1.00 904 1824
## wgPMG -0.26 0.09 -0.43 -0.08 1.00 910 1813
## wgPMM -0.07 0.09 -0.24 0.11 1.00 907 1734
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.23 0.01 0.21 0.25 1.00 7758 10394
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
We will now look at the model’s predicts for the average participant (i.e, ignoring the random intercepts).
Fit Model to Psychophysical Data
We will now fit the model to the data.
## Family: lognormal
## Links: mu = identity; sigma = identity
## Formula: threshold ~ wg - 1 + (1 | subject)
## Data: d_dispthresh (Number of observations: 186)
## Samples: 4 chains, each with iter = 10000; warmup = 5000; thin = 1;
## total post-warmup samples = 20000
##
## Group-Level Effects:
## ~subject (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.39 0.11 0.24 0.66 1.00 2870 4657
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## wgCM -0.22 0.16 -0.53 0.11 1.00 2403 4638
## wgCMM -1.15 0.16 -1.47 -0.82 1.00 2306 4800
## wgP2 -0.29 0.17 -0.62 0.06 1.00 2390 4846
## wgP3 -0.69 0.16 -1.00 -0.36 1.00 2419 4513
## wgP31M -1.18 0.16 -1.49 -0.85 1.00 2299 4123
## wgP3M1 -1.34 0.16 -1.66 -1.02 1.00 2360 4265
## wgP4 -0.89 0.16 -1.21 -0.56 1.00 2458 4682
## wgP4G -1.22 0.17 -1.54 -0.88 1.00 2417 4194
## wgP4M -1.29 0.16 -1.60 -0.96 1.00 2305 4474
## wgP6 -1.20 0.17 -1.52 -0.86 1.00 2489 4469
## wgP6M -1.41 0.17 -1.74 -1.07 1.00 2248 5050
## wgPG 0.36 0.17 0.04 0.71 1.00 2305 4597
## wgPGG -0.31 0.17 -0.62 0.02 1.00 2327 4833
## wgPM -0.79 0.17 -1.11 -0.44 1.00 2457 4343
## wgPMG -0.96 0.16 -1.28 -0.64 1.00 2330 4618
## wgPMM -1.17 0.16 -1.49 -0.84 1.00 2355 4413
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.42 0.02 0.37 0.46 1.00 10723 13316
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
EEG Control Data
We will also fit models to the control data. As we can see from Figure 2.4, the group differences are much smaller.
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#bulk-ess